Introduction

GravitationalWaves: a Python package designed to simulate, detect, and analyze continuous gravitational wave-forms. In addition to creating simulations of gravitational waves, the package also takes in observed data for comparison in detection and analysis.

GravitationalWaves is a Python package for simulating the gravitational wave-forms of binary mergers of black holes and neutron stars, computing several properties of these mergers and wave-forms, and evaluating their detectability. In addition, the package also takes in observed data from gravitational wave detectors, such as the Laser Interferometer Space Antenna (LISA), to compare data and predict detection rates for stellar-origin binaries.

The Laser Interferometer Space Antenna (LISA) is a proposed space probe to detect and accurately measure gravitational waves—tiny ripples in the fabric of spacetime—from astronomical sources. LISA would be the first dedicated space-based gravitational wave detector.

This report details software design decisions and analysis of the GravitationalWaves package system, with respect to existing software libraries. The software libraries used for the basis of this project include legwork and riroriro. Legwork is a python package for computing evolution and detection rates of gravitational-wave sources discovered by space-based detectors such as LISA. Riroriro is a python package for simulating and evaluating gravitational waves.

In the following, I address the major software design decisions, software requirements, usage guidelines, and the extensibility of the software for the GravitationalWaves package.


Background

Gravitational waves, first theorized by Einstein, arise during events such as the collision of supermassive binary black hole mergers. A binary black hole is a system consisting of two black holes in close orbit around each other. These collisions are so powerful that they create distortions in space time, known as gravitational waves.

Definition. Gravitational waves are invisible distortions in space time, caused by massive events such as collisions between two black holes or neutron stars.

The Laser Interferometer Space Antenna is a proposed space probe to detect and accurately measure gravitational waves—tiny ripples in the fabric of spacetime—from astronomical sources. LISA would be the first dedicated space-based gravitational wave detector.

The space-based gravitational wave detector LISA presents a new view of gravitational waves by focusing on lower frequencies \(\left( 10^{−5} \lt f / \mathrm{Hz} \lt 10^{-1} \right)\) than ground-based detectors. This enables the study of many new source classes including mergers of supermassive black holes, extreme mass ratio inspirals, and cosmological GW backgrounds.

The importance of studying gravitational wave-forms stems from the idea of detecting and using gravity to estimate other dynamical astrophysical phenomena, giving enormous potential for discovering parts of the universe that are invisible to the eye, such as black holes and other unknowns.


The Package

GravitationalWaves is a Python package for simulating the gravitational wave-forms of binary mergers of black holes and neutron stars, computing several properties of these mergers and wave-forms, and evaluating their detectability. In addition, the package also takes in observed data from gravitational wave detectors, such as the Laser Interferometer Space Antenna (LISA), to compare data and predict detection rates for stellar-origin binaries.

Specifically, we implement the expressions for the evolution of binary orbits due to the emission of gravitational waves, equations for the strain amplitudes and signal-to-noise ratios of binaries, and approximations for the LISA sensitivity curves.

Package Data

GravityPy uses simulated data for gravitational waves and observations from detectors such as the Laser Interferometer Space Antenna (LISA).

Simulated Data: We obtain the simulated data by implementing the package functions for constructing simulations of gravitational waveform signals from a binary merger of two black holes or two neutron stars and outputting the data of such signals in terms of frequency and strain amplitude.

i_amp = ins.inspiral_strain_amplitude(i_Aorth, i_Adiag)
i_freq = mat.frequency_SI_units(i_m_omega, M)
Strain Amplitude Frequency
7.9097e-22 10.0066
7.9105e-22 10.0081
7.9112e-22 10.0096
7.9121e-22 10.0111
7.9128e-22 10.0126
7.9136e-22 10.0141
7.9144e-22 10.0156

Observed Data: The observed data contains LISA-validated binaries obtained from the Kupfer et al. 2018 study. These observations are preloaded to the package and used to help generate the Source class based on these provided binaries.

Source f (mHz) \(\beta\) (mas) \(\sigma_\beta\) (mas) d (pc) \(\sigma_d\) (pc) \(\mathcal{A}\) SNR
AM CVn type
HM Cnc 6.22 NA NA [5000] 6.4 211.1 ± 3.18
V407 Vul 3.51 0.095 0.327 1786 667 11.0 ± 5.9 169.7 ± 2.17
ES Cet 3.22 0.596 0.108 1584 291 10.7 ± 4.6 154.3 ± 2.09
SDSSJ135154.46-064309.0 2.12 0.596 0.313 1317 531 6.2 ± 3.5 21.8 ± 0.24
AMCVn 1.94 3.351 0.045 299 4 28.3 ± 3.2 101.2 ± 0.96
SDSS J190817.07+394036.4 1.84 0.954 0.046 1044 51 6.1 ± 2.4 20.3 ± 0.13
HP Lib 1.81 3.622 0.052 276 4 17.5 ± 3.9 43.7 ± 0.28
PTFl J191905.19+481506.2 1.48 0.550 0.327 1338 555 3.2 ± 1.8 4.0 ± 0.02
CXOGBS Jl 75107.6-294037 1.45 1.016 0.146 971 156 4.2 ± 1.8 4.5 ± 0.02
CR Boo 1.36 NA NA 337^a ±44^a-35 13.4 ± 4.2 21.9 ± 0.13
V803 Cen 1.25 NA NA 347^a ±32^a-27 16.0 ± 5.4 26.2 ± 0.17
Detached white dwarfs
SDSS J065133.34+284423.4 2.61 1.000 0.476 933 493 16.2 ± 8.6 90.1 ± 1.13
SDSS J093506.92+441107.0 1.68 NA NA 645^b 41^b 29.9 ± 7.7 44.9 ± 0.31
SDSS J163030.58+423305.7 0.84 0.937 0.270 1019 357 11.5 ± 4.9 4.6 ± 0.03
SDSS J092345.59+ 302805.0 0.51 3.340 0.173 299 10 26.4 ± 6.5 5.6 ± 0.06
Hot subdwarf binaries
CD-30°11223 0.47 2.963 0.080 337 9 41.5 ± 1.8 4.9 ± 0.04

Limitations: Both simulated and observed calculations apply a low-order post-Newtonian description of gravitational wave emission, meaning that the data do not account for higher-order effects.

Software Requirements

programming language

Package Version Type Description
python \(\geq\) 3.7.0 Programming language High-level, general-purpose programming language
pip \(\geq\) 21.0.0 Python package Recommended package installer for Python
numba \(\geq\) 0.5.0 numba NumPy-aware optimizing compiler using LLVM
numpy \(\geq\) 1.16.0 numpy Fundamental package for array computing
astropy \(\geq\) 4.0.0 astropy Astronomy and astrophysics core library
scipy \(\geq\) 1.5.0 scipy Software for mathematics, science, and engineering
matplotlib \(\geq\) 3.3.2 matplotlib Library for creating various visualizations in Python
seaborn \(\geq\) 0.11.1 seaborn Library for making statistical graphics in Python
schwimm \(\geq\) 0.3.2 schwimm Common interface for parallel processing pools
legwork \(\geq\) 0.2.4 legwork Package for analyzing gravitational waves

List of Dependencies:

Project Structure

The structure of the GravitationalWaves github repository is as follows:

gravitational-waves
 ┣ docs (folder contains functional and design specification)
 ┃ ┗ presentation.Rmd
 ┣ examples
 ┃ ┣ 01_InstantiateSourceClass.ipynb
 ┃ ┣ 02_CalculateSNR.ipynb
 ┃ ┣ 03_PlotSourceDistribution.ipynb
 ┃ ┣ 04_Visualizations.ipynb
 ┃ ┗ 05_SimulateGravitationalWaves.ipynb
 ┣ GravitationalWaves (python package folder)
 ┃ ┣ tests
 ┃ ┃ ┣ test_psd.py
 ┃ ┃ ┣ test_source.py
 ┃ ┃ ┣ test_strain.py
 ┃ ┃ ┣ test_utils.py
 ┃ ┃ ┣ test_visualization.py
 ┃ ┃ ┣ test_wavesim.py
 ┃ ┃ ┗ __init__.py
 ┃ ┣ psd.py
 ┃ ┣ source.py
 ┃ ┣ strain.py
 ┃ ┣ utils.py
 ┃ ┣ visualization.py
 ┃ ┣ wavesim.py
 ┃ ┗ __init__.py
 ┣ LICENSE (open source MIT License)
 ┣ README.md (overview of the project)
 ┗ setup.py (file initializes the project after it has been cloned)

From the above, we can see the top-level package called GravitationalWaves, with six modules: strain.py, psd.py, utils.py, visualization.py, source.py, wavesim.py. There’s also the subpackage called tests located in the top-level package. There’s also an examples and docs directory that contains demos and documentation of the package, though, they’re not packages since both don’t’t include an __init__.py.


Software Components

The software components for the GravitationalWaves package include the package’s modules, classes, and functions. A Python package is any directory with an __init__.py file, and this file gathers all package-wide definitions to import every module in the package. Python modules are an essential abstraction layer, which allow separating code into parts holding related data and functionality. For example, a layer of a project might handle interfacing with user actions, while another handles data manipulation. To separate these two layers, we’d regroup all user interfacing functionality into one file and all data manipulation functionality into another file.

Detailed specification (ComponentSpec)

For the following component specification, we define each components’ metadata. interface, and implementation. The metadata consists of a component’s name and description. The interface specifies the list of inputs/outputs and their properties. And the implementation specifies how to execute the component instance, pass inputs to the component code, and get the component’s outputs. After the component specification, I introduce the component code, which implements the logic needed to perform a step in the package workflow.

Module 1. strain.py

Name: strain

URL: http://www.pobox.com/~tranter

Description: A program that uses functions from LEGWORK to calculate signal-to-noise ratio of sources based on four different cases. For the different cases, the module contains four functions that cover whether a source is circular-stationary, circular-evolving, eccentric-stationary, or eccentric-evolving.

Metadata: Standard object’s metadata:

  • Annotations: A string key-value map used to add information about the component. Currently, the annotations get translated to Kubernetes annotations when the component task is executed on Kubernetes.

inputs: - {name: Training data name of the input/output}

  • {name: Rounds, description: ‘of the input/output’, default: ‘30’, type: (String, Integer, Float, and Bool), optional: Specifies if input is optional or not (TRUE FALSE)}

outputs: - {name: Trained model, type: XGBoost model, description: ‘Trained XGBoost model’}

implementation: container: image: :b3a64d57 command: [ /ml/train.py, –train-set, {inputPath: Training data}, –rounds, {inputValue: Rounds}, –out-model, {outputPath: Trained model}, ] ““” Computes strain amplitude at n harmonics

### Parameters:
    m_c (`float/array`): Chirp mass of each binary.
    f_orb (`float/array`): Orbital frequency of each binary at each timestep.
    ecc (`float/array`): Eccentricity of each binary at each timestep. 
    n (`int/array`): Harmonic(s) at which to calculate the strain. 
    dist (`float/array`): Distance to each binary. 
    position (`SkyCoord/array`, optional): Sky position of source. Defaults to None.
    polarisation (`float/array`, optional): GW polarisation angle of the source. Defaults to None.
    inclination (`float/array`, optional): Inclination of the source. Defaults to None.
    interpolated_g (`function`, optional): Computes g(n,e) from Peters (1964). Defaults to None.

### Returns:
    h_0 (`float/array`): Strain amplitude. Shape is (x, y, z).
"""

Module 2. psd.py

Name: psd

URL: http://www.pobox.com/~tranter

Description: The psd module evaluates the effective noise power spectral density of a detector at different frequencies. The module contains the LISA sensitivity curve that can be tweaked by adjusting parameters such as the observation time, response function, and the arm length. For the Galactic confusion noise, the module uses the Robert et al. 2019 model.

Metadata:

  • Annotations:

inputs: - {name: Training data name of the input/output}

  • {name: Rounds, description: ‘of the input/output’, default: ‘30’, type: (String, Integer, Float, and Bool), optional: Specifies if input is optional or not (TRUE FALSE)}

outputs: - {name: Trained model, type: XGBoost model, description: ‘Trained XGBoost model’}

implementation:


Module 3. utils.py

Name: utils

URL: http://www.pobox.com/~tranter

Description: The utils module is a collection of miscellaneous utility functions for computing conversions between variables, the chirp mass of binaries, and expressions from the Peters equation (Peters and Mathews 1963), such as the relative power of gravitational radiation and enhancement factor.

Metadata: Standard object’s metadata:

  • Annotations:

inputs: - {name: Training data name of the input/output}

  • {name: Rounds, description: ‘of the input/output’, default: ‘30’, type: (String, Integer, Float, and Bool), optional: Specifies if input is optional or not (TRUE FALSE)}

outputs: - {name: Trained model, type: XGBoost model, description: ‘Trained XGBoost model’}

implementation:


Module 4. visualizations.py

Name: visualizations

URL: http://www.pobox.com/~tranter

Description: The visualization module uses matplotlib and seaborn to plot quick visuals for analyzing of a collection of sources. The plots include 1- and 2-dimensional distributions and the LISA sensitivity curve.

Metadata:

  • Annotations:

inputs: - {name: Training data name of the input/output}

  • {name: Rounds, description: ‘of the input/output’, default: ‘30’, type: (String, Integer, Float, and Bool), optional: Specifies if input is optional or not (TRUE FALSE)}

outputs: - {name: Trained model, type: XGBoost model, description: ‘Trained XGBoost model’}


Module 5. source.py

Name: Source (class)

URL: http://www.pobox.com/~tranter

Description: The source module is the central module of the package, which provides a simple interface to the functions in the remaining modules. This module contains the Source class for analyzing a generic set of gravitational wave sources, stationary/evolving and circular/eccentric.

Metadata: Standard object’s metadata:

  • Annotations:

inputs: - {name: m_1, description: , type: Integer, optional: FALSE} - {name: m_2, description: , type: Integer, optional: FALSE} - {name: ecc, description: , type: Integer, optional: FALSE} - {name: dist, description: , type: Integer, optional: FALSE} - {name: n_proc, description: , type: Integer, default: ‘1’, optional: TRUE} - {name: f_orb, description: , type: Integer, optional: FALSE} - {name: a, description: , type: Integer, default: None, optional: TRUE} - {name: position, description: , type: Integer, default: None, optional: TRUE} - {name: weights, description: , type: Integer, default: None, optional: TRUE} - {name: gw_lum_tol, description: , type: Integer, default: 0.05, optional: TRUE} - {name: stat_tol_e, description: , type: Integer, default: 1e-2, optional: TRUE} - {name: interpolate_g, description: , type: Integer, default: True, optional: TRUE} - {name: interpolate_sc, description: , type: Integer, default: True, optional: TRUE}

outputs: - {name: Trained model, type: XGBoost model, description: ‘Trained XGBoost model’}

implementation:

import GravitationalWaves
import GravitationalWaves.source 
import astropy.units as u

source1 = source.Source(m_1 = 11 * u.Msun, m_2 = 11 * u.Msun, ecc = 0.3, 
                        dist = 9 * u.kpc, f_orb = 1e-4 * u.Hz, interpolate_g = False)

Module 6. wavesim.py

Name: wavesim

URL: http://www.pobox.com/~tranter

Description: The gwsim module contains functions for simulating gravitational waves from binary black holes. The procedure includes parts for simulating the inspiral portions, pre-matching and post-matching parts for simulating the merger/ringdown portions, and parts for matching together the inspiral and merger/ringdown portions.

Metadata: Standard object’s metadata:

  • Annotations:

inputs: - {name: Training data name of the input/output}

  • {name: Rounds, description: ‘of the input/output’, default: ‘30’, type: (String, Integer, Float, and Bool), optional: Specifies if input is optional or not (TRUE FALSE)}

outputs: - {name: Trained model, type: XGBoost model, description: ‘Trained XGBoost model’}

implementation:


Use Cases (functional spec)

With this package, users can

  • Simulate a gravitational wave signal from a binary merger of two black holes, and output the frequency and strain amplitude data for that signal.
  • Create a random collection of possible LISA binary sources and calculate their strain amplitudes given any range of frequency harmonics.
  • Calculate and plot the sensitivity curve of a gravitational wave detector
  • Check whether a binary source is eccentric/circular and evolving/stationary
  • Compute the signal-to-noise ratio (SNR) of a single source binary system, as well as a collection of systems
  • Plot either 1D or 2D parameter distributions for a collection of sources
  • Determine the number of detectable binaries from a population of sources

 

Demo 1. Single source SNR calculation

The most basic use case for GravitationalWaves is to calculate the signal-to-noise ratio of a single stellar-mass binary system. We do this by using the package’s source module to generate a Source class. First we must import the package module. As soon as we use import statements, we can use modules, which can be built-in, installed third-party modules, or internal project modules. The best import statement is simply import modu. Using from modu import func is a way to pinpoint the function to import and put it in the local namespace.

To import the file source.py in a directory, we use the statement import GravitationalWaves.source. This statement looks in the directory for __init__.py and then for source.py, executing all top-level statements in both files. Consequently, any variable, function, or class defined in source.py is available in the local namespace.

import GravitationalWaves
import GravitationalWaves.source 
import astropy.units as u

source1 = source.Source(m_1 = 11 * u.Msun,
                        m_2 = 11 * u.Msun,
                        ecc = 0.3,
                        f_orb = 1e-4 * u.Hz,
                        dist = 9 * u.kpc,
                        interpolate_g = False)
                          
source1.get_snr()

For this example, GravitationalWaves checks whether the source is eccentric/circular and evolving/stationary, and chooses the fastest method to accurately calculate the SNR.

Demo 2. Multiple source SNR calculate

In the next example, we use GravitationalWaves to calculate the detectability of a collection of sources. We first import the package, then create a random set of possible LISA sources.

import GravitationalWaves.source as source
import numpy as np
import astropy.units as u

# create a random collection of sources
nVals = 1800
mass1 = np.random.uniform(0, 12, n_values) * u.Msun
mass2 = np.random.uniform(0, 12, n_values) * u.Msun
eccInit = 1 - np.random.power(6, n_values)
lumDist = np.random.normal(9, 1.5, n_values) * u.kpc
orbFreq = 12**(-6 * np.random.power(3, n_values)) * u.Hz

Using these random sources, we can instantiate a Source class.

sources = source.Source(m_1 = mass1, m_2 = mass2, 
                        ecc = eccInit, dist = lumDist, 
                        f_orb = orbFreq)

Now, we can calculate the SNR (signal-to-noise ratio) for these sources. This function splits the sources based on whether they are stationary/evolving and circular/eccentric.

snr = sources.get_snr(verbose = True)

The SNR values are now stored in sources.snr, and we can mask any values that don’t meet some detectable threshold. In the following, we set the threshold to 7.

detectableSources = sources.snr > 7
print("{} out of the {} sources are detectable".format(
  len(sources.snr[detectableSources]), nVals))

Comparison

Comparison to the design of some existing software library


Discussion

Software Design Decisions:

Overall, this project taught me how to structure a Python project. Structuring GravityPy means making decisions about how the project can best achieve its goals: simulating, analyzing, and detecting gravitational waves. In other words, the structure requires writing clean code, with clear logic and dependencies, and organizing files and folders in the file system. Writing clean code with clear logic is part of the architectural task of building out the different components of a project and their interactions. Organizing a project requires making decisions, such as which features should go into a module, how data should flow through the project, and which features and functionality should be grouped together.

By structuring the Python project, I learned how to divide the code base into clean, efficient modules using packages. This step requires building modules as an abstraction layer to divide the code into sections containing related data and functionality. For example, this project separates the layers that deal with plotting results from the layers that deal with miscellaneous utility functions by grouping all plotting functions in the visualization module and all utility functions in the utils module. Overall, the straightforward module importing model and packaging system provided in Python makes it notably easier to structure a Python package.

Discussion of at least two design decisions made and why you made them 6

Design considerations: * Giving optional parameters for more complex stuff by user. Extensibility! * Forces users to be using Pandas DF.

scikit-learn workflow: from sklearn.ensemble import RandomForestClassifier from sklearn.datasets import load_iris from sklearn.preprocessing import test_train_split

data = load_iris() model = RandomForestClassifier(n_trees=10) train_X, train_y, test_X, test_y = test_train_split(data.drop([“label”], axis=1), data[“label”]) model.fit(data)

  • Maintain information about the transformations that take place.

Software Extensibility.

Discussion of extensibility of your software (i.e., how difficult is it to extend?)

Overall, GravityPy is designed to help calculate binary sources of gravitational waves, whether simulated or LISA-like. This package aims to provide a way to study and better understand the detectability of such compact object binaries. Future work includes adding more functions, equations, and modules to the package to implement and analyze gravitational-wave emission, gravitational wave strain, SNR, and visualization modules, and see these implementations’ effect on orbital evolution. This will require adding more modules and functions to compute higher-level linear algebra-based equations and mathematical models.


References

Ashton, Gregory, David Keitel, Reinhard Prix, and Rodrigo Tenorio. 2020. PyFstat.” Zenodo. https://doi.org/10.5281/zenodo.3967045.
van Zeist, Wouter G. J., Héloïse F. Stevance, and J. J. Eldridge. 2021. Riroriro: Simulating gravitational waves and evaluating their detectability in Python.” The Journal of Open Source Software 6 (59): 2968. https://doi.org/10.21105/joss.02968.
Wagg, Tom, Katelyn Breivik, and Selma E. de Mink. 2021. LEGWORK: A python package for computing the evolution and detectability of stellar-origin gravitational-wave sources with space-based detectors.” arXiv e-Prints, November, arXiv:2111.08717. https://arxiv.org/abs/2111.08717.